#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <time.h>
#include <stdio.h>
#include "ran_tg.h"
#define NSIMS 10*1000*1000

double ran_tg_robert(gsl_rng *r,double a)
{
    double t,logu,z;

    if (a<0.)
        return ran_tg_reject(r,a);

    t = 0.5*(a+sqrt(a*a+4.));
    do
    {
        z = a - log(gsl_rng_uniform(r))/t;
        logu = log(gsl_rng_uniform(r));
    }
    while (logu>-0.5*gsl_pow_2(z-t));
    return z;
}

double ran_tg_devroye_reject(gsl_rng *r,double a)
{
    if (a<0.65)
        return ran_tg_reject(r,a);
    else
        return ran_tg_devroye(r,a);
}

double ran_tg_inverse(gsl_rng *r,double a)
{
    double u,z;

    if (a<5.) {
      u = gsl_rng_uniform(r);
      z = gsl_cdf_ugaussian_P(-a);
      return -gsl_cdf_ugaussian_Pinv(z*u);
    }
    else
      return ran_tg_devroye(r,a);
}


